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Abstract 

This paper introduces a method to preform optical tomography, using 3D ra¬ 
diative transfer as the forward model. We use an iterative approach predicated on 
the Spherical Harmonics Discrete Ordinates Method (SHDOM) to solve the opti¬ 
mization problem in a scalable manner. We illustrate with an application in remote 
sensing of a cloudy atmosphere. 


1 Introduction 

Optical tomography is a 3D imaging technique that uses optical measurements on the 
boundary of a domain to find the spatial distribution of parameters within. It has nu¬ 
merous applications in bio-medical imaging. For a list of applications we refer the 
reader to EH) and references therein. In sharp contrast, operational remote sensing of 
clouds in the Earth’s atmosphere using space-borne optical sensors is predicated on the 
classic plane-parallel medium geometry with horizontally uniform properties, hence 
the ID radiative transfer equation (RTF); this reduces the problem to estimating a sin¬ 
gle cloud optical thickness, and shifts the focus toward microphysical properties ffl- 
Some of the biases in retrieved cloud properties caused by this sometimes very crude 
approximation of real cloud structure are documented in 

Solving the inverse problem using the full 3D RTF as a forward model can be dif¬ 
ficult and computationally demanding. For some applications, it is possible to use an 
approximate model. In dense media (mean free path small compared the overall size), 
with scattering dominating absorption, it is possible to use a diffusion approximation 
of the 3D RTF. This results in the inverse problem of Dijfuse Optical Tomography 
(DOT) |[7}|^. When the mean free path is large compared to the outer size of the 
medium, the measured radiant energy is dominated by direct and single scattered in¬ 
tensities. The resulting inverse problem is single scattering tomography |[2l p^pH| . 
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Figure 1: (a) Domain and boundary conditions; (b) Aperture function. Blue marks the 
spatial support and red marks the angular support 

Other approximations and their derivations can be found in (13 We propose to solve 
the inverse problem using the RTF, with neither single scattering, nor diffusion ap¬ 
proximations. However, to make the inverse problem tractable, we derive an iterative 
optimization framework that uses different parts of the Spherical Harmonics Discrete 
Ordinates Method (SHDOM) 3D RTF solver p3| . 

2 Theoretical background 

2.1 Radiative transfer 

Our forward model is steady-state 3D radiative transfer. This model is used in passive 
imaging when source gating is sufficiently slow. The 3D RTF governs the transport 
of light through a medium with distribution of scattering and/or absorption. Consider 
a domain ft having boundary dft whose outward facing normal is (Fig. [^). The 
transport domain is indicated by position x G and direction of propagation u? G 
(unit sphere). The radiation (light) field is Ix The subscript A indicates 

wavelength dependency and will be omitted from here on, for simplicity. The boundary 
condition (Fig.[^) is 


/ (CC, LV) = /incident {x,(jj) if (jJ ■ 'd < 0, CC G Oft, 


( 1 ) 


where u; • < 0 defines incoming radiation, being the outgoing normal to dft. The 

transport of light is formulated in terms of the following conservation law p^p~5| : 


u) • V/ (cc, (jij) = /3 (x) [J (tc, (jj) — I (cc, cj)] X E ft, 


( 2 ) 


where /3 is the extinction coefficient. Here J {x,uj) is the source function. In the 
absence of emission within the medium, it is also known as the in-scattering term. 
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since it accounts of an increase of radiation due to in-scattering. In this case, it is 
expressed as 


J (cc, u) 


w 

47r 



u) • uj') I (cc, u)') du?', 


( 3 ) 


where w is the single scattering albedo andp (u; • Cc;'|^is the phase function. The phase 
function describes the probability of a photon traveling in direction u' into scatter to 
direction u). The phase function satisfies a normalization condition 
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u) • uj') du)' = 1. 


( 4 ) 


Eqs. ([^-([^, and the boundary condition Q, define a complete radiative transfer for¬ 
ward model for an externally illuminated medium. 

Integrating Eq. ^ along a specific direction lv results in the integral form of the 
3D RTE |[T4l[T^ 




Jo exp 



dr 


+ 


Xq 


j J{x',u})l3 {x') 


exp 



dx\ 


( 5 ) 


bearing in mind the expression of J {x^uj) in terms of / {x^uj) in Here Xq is 
a point on the boundary and Iq holds the boundary condition Q- The operation we 
denote as 


j f{r)dr (6) 

X 

is simply a ID (line) integral over a field f{x) along the segment extending from x to 
x'. Numerically, this is preformed by hack-projecting (BP) a ray through the medium. 

The expression for / (cc, cj) in for a given J (cc, u;), is often referred to as the 
up-wind sweep of radiant energy sources. We adopted the Spherical Harmonics Dis¬ 
crete Ordinates Method (SHDOM) (MI) for the forward 3D RT, which is available as 
open-source code. In SHDOM, ^ is indeed used in a straightforward post-processing 
stage for estimating the observed quantity I (x^lv) after a nontrivial numerical pro¬ 
cedure is performed to obtain J {x^uj) throughout the grid. This staged approach is 
exploited further on. We now formalize our forward model in operator language. 


2.2 Operator Notation 

We follow the definitions of to formulate the forward model using operator nota¬ 
tions (Table summarizes the notations). Denote by 0 the space of extinction fields 

4n principle, the phase function depends on both incoming and outgoing directions (u;, u;'). However, 
dependency is often assumed to be solely on the scattering angle, equivalent to u; • u;'. 
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Operator 

Symbol 

Domain 

Range 

Radiance forward mapping 

X 

0 - Extinction fields 

Z - Radiance fields 

Measurement 

M 

Z - Radiance fields 

y - Measurements 

Forward operator 

T 

0 - Extinction fields 

y - Measurements 

Scatter forward mapping 

J 

0 - Extinction fields 

V - Scatter fields 

Transformation operator 

T 

V - Scatter fields 

Z - Radiance fields 


Table 1: Operator notation summary 


over the domain ft. Let ^ (x) G 0 be an extinction field over the domain ft. Let S 
represent a set of radiation sources over dfl. For a given source s e S, denote Zs as the 
set of all possible radiation fields that satisfy Eqs. Ql' across all possible /3 e &. The 
set Zs is infinite, since for each extinction field P{x) there is a corresponding radiation 
field Is = Zs (P). More generally the radiance forward mapping, X : 0 ^ is a 
mapping from the optical parameter domain, to the radiance field range Z = 
Measurements are an operator Ai : Z ^ y, mapping the continuous function space 
to a vector space y. Consequently the forward operator X" is defined as 

X = MX lO^y. (7) 

The measurement operator M is defined by the detector’s aperture function, w e Q x 
(Fig#). The aperture function defines the manner in which the detector collects radi¬ 
ance, over a spatial and angular support. A given aperture and source pair (w,s) yields 
an element 


yw,s = ^W,S (/?) = Myjis (/ 3 ) = {w, Is)n, ( 8 ) 

for a specific extinction field P, where 

(•, •)n = j J ■ ■ dudx. (9) 

For an idealized single-pixel detector positioned at cc*, collecting radiation flowing in 
direction u;*. 


yw,s = {6{x- X*) (5 (u; - u;*) , = h . (10) 

We can define a matrix Y whose elements y^^^s correspond to different source-detector 
configurations. A column of Y represents measurements by a single detector, for mul¬ 
tiple sources. A row represents measurements by multiple detectors, for one particular 
light source. The vector y is the column stack of Y . 


2.3 Optical tomography 


Using the operators defined in Sec. 2.2 we express the process of optical tomography. 
Tomographic reconstruction is an estimator of /3 that minimizes a defined cost 


/3 = arg min {£ [y, T (/3)] + (/3)} . 


( 11 ) 


4 






0(x) 



b{x) 



Figure 2: Illustration of a ID linear interpolation kernel 


where £ {/3)] is the data fit (fidelity) functional and ^ (;5) is a regularization on 

the optical parameters. Here a is a tunable parameter, chosen in accordance with the 
noise level, to balance the two terms. 


Except for a few special geometric configurations mg. 3D tomographic recovery 
cannot be done analytically. Thus, the continuous function is often discretized (Fig.|^ 

ATgrid 

P{x)= ^ Pkbk (x), (12) 

k=l 

where {f3k}k=i^ discrete parameters, bk {x) is an interpolation kernel and A^grid is 
the number of grid points used for discretization. We apply ( p^ to Eq. ( pT] ) to seek an 
estimator for (5 = (/^i,..., , where (• denotes transposition. 


We focus here on gradient-based optimization methods. The forward operator T 
is a non-linear function of the optical parameters. Thus, one approach Q estimates 
(11) by linearizing T. Suppose the solution is = P^ ^ 5P, which is a perturbation 
of an initial guess P^. It is possible to linearize T and solve dp using a linear set for 
equations. However, this approach requires an initial guess very close to the true so¬ 
lution. Another approach fT9}j2^ iteratively estimates the gradient with respect to (3. 
However, this approach has high computational complexity, as it requires O (A/'grid) 
simulations of the forward model per iteration (i.e., for a single computation of the 
gradient). Forward model simulation is time consuming particularly in the presence of 
multiple scattering. Simulating the forward model per element of the gradient does not 
scale well as A'grid increases. 


Our approach is also iterative, however each iteration is computationally simple, 
having run-time independent of A'grid- The approach does not directly optimize the 
nonlinear T. Instead, X is decomposed into a product of two operators. Each is opti¬ 
mized in time. This concept is analogous to expectation minimization (EM) optimiza¬ 
tion, which is used in DOT | |23l - [^ . 
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Figure 3: (a) A line integral over the in-scatter field (b) Integration over the field / 
for an arbitrary aperture spatial and angular support w 


3 Optimization Approach 

3.1 Forward map decomposition 

In Sec. 1^ the forward operator was defined in terms of the radiance field I. We now 
show that defining T in terms of the in-scatter field J considerably speeds up the 
computation time. We decompose the radiance forward mapping into two operators, 
T = TJ, which we now introduce. 

For a particular radiation source 5 G 5, let Js = be the in-scatter field that 

satisfies Eq. while Is satisfies Eqs. Define as the set of possible in-scatter 

fields for a particular source. The in-scatter forward mapping, : 0 ^ V, is a 
mapping from the optical parameter domain to the in-scatter field range V = Us 
We define the following transformation operator 

Is=T{p)Js. (13) 

Eq. ^ defines a relation between a given light field / and a corresponding in-scatter 
field J. The transformation operator 'T, defined by Eq. defines the inverse relation: 
for a given in-scatter field J, a corresponding light field / is attained. Consequently, 
we define the forward operator in terms of the in-scatter forward mapping 

T = Mr{p)j{p):e^y. (U) 


Define ^ as a line of sight from x* in direction —u;* (Fig.[^). For a convex spatial do¬ 
main, the intersection point between line of sight and the domain’s boundary is unique 
and denoted hy Xq = dfl D i (Fig.l^). 

Operator T is defined by Eq. For source s G 5, the response of an idealized 
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Figure 4: (a) Gradient decent optimization, (b) Surrogate function iterative optimiza¬ 
tion. Both depend on the initial guess Pq 

single pixel with aperture w = ^(cc — — u;*)is 

Vw^S — Js 



~ Xo 

n 

Xq 

n 

~ X 

n 

Js(a;o,w*)exp 

— 13 (r) dr 

J 

J 

— 13 {r) dr 

J 


- X* 

X* 

- X* 


( 15 ) 

Eq. ([Tg is simply an accumulation of the scattered radiance along the line of sight 
weighted by its corresponding attenuation factor (Fig. |^). Applying a general mea¬ 
surement operation is a manner of integrating over the aperture function of the detector 
(Fig#). 

3.2 Iterative surrogate function optimization 

Optimizing Eq. ^TT\ directly over p is computationally expensive. Finding the gradient 
direction in each iteration requires O (A/’grid) numerical computations of the forward 
model. With our approach, however, we iteratively keep Jg constant and estimate (3 
solely based on 'T. Instead of optimizing through a function that is difficult to compute 
(Fig.|^) we optimize /3 using a surrogate function p6| , which is efficiently computed 
(Fig.#). We define the following iterative optimization process. Define Pn as an 
estimate of p in the “n” iteration. The first step consists of computing the in-scatter 
field, which corresponds to the current estimate ( 3 n 

Jn= J iPn) . 

In the second step, keeping the in-scatter field constant, we solve the following opti¬ 
mization problem to find our next estimate of p 

/3„+i =argnnn{£:[j/,:F„(/3)]+a4'(/3)}, with (16) 

£: [y, (/?)] = [y- MTiP) Jnf [y - Jn ], 
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Figure 5: (a) Each gradient element is computed by tracing all the rays. Only the red 
rays contribute to the computation of the red grid point gradient, (b) A single trace of 
each ray is performed. While tracing each ray we compute the gradient contribution to 
each relevant grid point 


where 


Tn{P)=Mr{P)Jn. (17) 

is the surrogate function, and I^meas is the covariance matrix of our measurements. 
In the remainder, we will not make use of regularization (set a = 0). 


Eq.([^ defines an iterative optimization process, 

i . Start with an initial guess /^q. 

ii . Based on the current estimate /3n, numerically find the in-scatter field Jn = J' {Pn)- 

i i i . Optimize (17) to find the next estimate 

i V. Repeat steps i i - i i i until convergence. 


3.3 Scalability 

For a given measurement vector y E y, wq solve the optimization of ( p^ with a 
gradient-based method. We use the discretization described in Assuming uncor¬ 
related measurements 

Smeas = diag (cr^eas)> ^Tmeas = Eh 

where A^meas = X Ns is the number of measurements, N^j being the number of 
detectors and Ns the number of sources. In this study based on radiometry, we use 
cFi = 0.03 X yi (si 3% noise due mostly to calibration error). 

The task is thus to recover the gridded extinction b’or this purpose we 

find the k'th gradient element. Without loss of generality we formulate the problem in 
terms of a single source and many detectors. 


d 

dpk 


£[y,:Fr,m 


Nn, 


E 

W = 1 


cr2 ^ “ 


. (/ 5 ) yw\ •Ni.y 




Jni 


(19) 
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For an ideal single-pixel detector positioned at cc*, collecting radiation flowing in di¬ 
rection u;* we get 


•A^y 




^w,k 


+ B, 


w,k 


( 20 ) 


where 


and 

^w,k 



~ Xq 

n 


~ Xq 

^w,k 

- bk (r) dr 

/(a;o, w*)exp 

— P (r) dr 

J 


J 


_ X* 

Xq 

r 

X 

X 


( 21 ) 


j J„{x,u}*) bk{x) - l3{x) j bk{r)dr exp -Jl3{r)dr 


dx. (22) 


Term A results from the boundary illumination and can be readily computed. Term 
B deflnes a line integral over a held, computed using back-projection of rays from 
the detector through the medium. A straightforward approach calculates each element 
of the gradient vector by summing all the back-projected rays from all the detectors 
(Alg.j^). The complexity of this approach is O (A/’meas x A/'grid) back-projecting op¬ 
erations. However, most of back-projected rays do not contribute to B since the inter¬ 
polation kernel bk{x) typically has a small support region (Fig.[^). Instead, for each 
point along a ray we deflne a support region of grid points, according to the support 
region of the interpolation kernel (Fig.|^). 

A speciflc ray will contribute to the computation of gradient elements associated 
with the grid points in the support region. In the case of a linear interpolation ker¬ 
nel, a support region is composed of the eight corner grid points that make a grid cell. 
Since the support region is typically small, for each point along the ray we only need to 
look at a flnite amount of neighboring grid points (Alg.[^). The complexity of this ap¬ 
proach is O (A/’meas) back-projccting operations. Hence, the computational complexity 
is independent of the number of grid points which makes this approach scalable. 


4 Application to remote sensing 

4.1 General setup 

We use large eddy simulation (LES) |27|28l to generate the microphysical quantities of 
a realistic cloud held (Fig.|^. We use Mie scattering theory (implemented in the avail¬ 
able code p^ ) to convert the microphysical quantities (liquid water content from LES, 
assumed effective droplet radius reff and effective variance of the droplet size distri¬ 
bution Vq^) into the required optical quantities (extinction coefficient, phase function, 
single scattering albedo). We take reff = 10/im and reff = 0.1 for Gamma-distributed 
droplet sizes, and the wavelength is set to 672 nm. Finally, we add the extinction and 
phase function due to Rayleigh scattering by air molecules; aerosols are ignored in this 
simulation. 
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1: for k = 1 ^ A^grid do 

2: Grad/c = 0 

3: for = 1 ^ A^meas dO 

4: BP to compute 


—“ ay 

5: Err = ^ [:F„ {(3) - y^] 

w 


meas 


3: Err = ^ [:F„ (P) - y^] 

4: for k G Support region do 

5: BP to compute 

6: Grad/c = Grad/^ + 


1: Grad = 0 

2: for IC = 1 ^ A^meas do 


meas 


6: Grad/c = Grad/^ + 


+ Err 


+ Err + By^^k) 


7: end for 

8 : end for 


7: end for 

8 : end for 


(a) Straightforward approach 


(b) Scalable approach 


Alg. 2: Comparison of the two approaches to compute the gradient of the surrogate 
function. Grad = (Gradi,GradATg^id) 

We seek to recover the extinction of the cloud droplets on a cartesian grid, where 
the air is taken to be a known parameter dependent only on the height. The boundary 
conditions for the domain are 

• Collimated solar radiation at the top of the atmosphere (TOA), at a zenith angle 


of 60°. 


• Open boundary for the side faces, to approximate an isolated cloud observed at a 
wavelength where the Rayleigh sky is almost black (Rayleigh optical thickness 
of the atmosphere is only 0.04532 at 672 nm). 

• Lambertian reflectance at the surface (Earth) with an albedo of 0.05, to mimic 
the dark ocean surface. 

At A = 672 nm, liquid water and air are non-absorbing, so the single scattering albedo 
w is everywhere unity. In order to And the in-scatter fleld we use SHDOM |T3| . 


4.2 SHDOM 


Running the forward model is a balance between speed and accuracy. Monte Carlo 
methods can handle very complex optical media. However, they compute radiomet¬ 
ric quantities by random sampling the inflnite domain of possible light paths. This 
introduces stochastic noise at the output, which can be controlled by increasing the 
number of samples (photon paths) and variance reduction techniques. When many ra¬ 
diometric outputs of the same scene are sought, as in the case of many viewpoints, a 
model that solves the RTE directly has a much more favorable runtime p2| . A dis¬ 
crete ordinates representation models the flow of radiant energy in the domain easily 
and intuitively The SHDOM model uses a grid for the spatial depen¬ 

dency, and a linear interpolation kernel bk{x). Spherical Harmonic expansion | p^[^ 
is used to compute angular integrals. SHDOM solves the forward model 3D RTE in 
an integral form for the source function J (cc, cj). It can adaptively truncate negligi¬ 
ble coefficients in the series expansion. It also has adaptive spatial grid capability and 
multi-grid acceleration. 
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Figure 6: Cumulus cloud field generated by the LE S This nadir-viewing image 
was generated with the MYSTIC Monte Carlo code |30||3H . We selected the isolated 
cloud circled in green to test our tomographic reconstruction algorithm. 


4.3 Simulation Results 

We simulate an atmosphere of size 2 km x 0.72 km x 1.44 km. The unknown extinc¬ 
tion is composed of 100 x 36 x 36 grid points (129,600 unknowns). The measurements 
are taken to approximate the resolution of AirMSPI (Airborne Multi-angle Spectro- 
Polarimeter Imager) p6| , which is an airborne sensor that can sample radiance at 10 m 
resolution and, to emulate the space-borne Multi-angle Imaging Spectro-Radiometer 
(MISR) at 9 viewing zenith angles: {±70.5°, ±60°, ±45.6°, ±26.1°, 0°} where 
± indicates two azimuthal half-planes at 180° apart (Fig. |7]). We simulate the mea¬ 
surements using SHDOM and add 3% gaussian noise to simulate the radiometric cali¬ 
bration noise, which is the dominant factor for this sensor. We initialize our algorithm 
with /3 = 0 (an atmosphere containing only air molecules). 

Simulated images as viewed from the AirMSPI instrument flying over the recov¬ 
ered cloud are shown in Fig. The converged reconstruction is displayed in Fig. 
along with the ground-truth. We compare the total extinction recovered, relative to the 
total extinction, with the LES/Mie-based truth. We define a relative error measure 

I ^ground truth _ ^recovered I 

Absolute Relative Error = ^ ^ ^ 

E ^ground truth 

k Pk 

The performance for the recovery is summarized in Table 

The extinction coefficient summed over the grid, Pk (in km“^), can be con¬ 
verted into a domain-average vertical optical thickness (x0.04 km/(36 x 100)), yield- 
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36 Grid points 


Figure 7: AirMSPI measurements: 9 zenith viewing angles (VZAs) in the along track 
direction. There are also rows of 36 across-track pixels (20 m width). 


ing 0.956. The recovered value is 1.089, a 14% overestimate. Most of this vertically- 
integrated opacity is of course in the cloud (the Rayleigh optical thickness of the lowest 
1.44 km layer being only 0.0075), which covers ^20% of the horizontal domain (based 
on LWC > 10“^g • m“^ threshold). 

5 Summary and Outlook 

We derive a novel iterative optimization method to perform radiative transfer tomog¬ 
raphy. The unknown /3 directly affects 'T{(3). From Eq. this principle can only 
retrieve parameters that relate to the extinction (we cannot estimate the phase func¬ 
tion using this formulation). Sec. |3.3| explains the computational advantage of this 
approach. We show an application to remote sensing of Earth’s atmosphere (Sec. [^, 
however, this approach could potentially be applied to preform optical tomography of 
biological tissues. We use SHDOM as our forward mapping engine. Nevertheless, it 
is possible to use different radiative transfer engines, such as Monte Carlo. Eurther 
enhancement of scalability may be obtained by use of adjoint operators | [38| , a sub¬ 
ject intended for further research. Another extension is to determine how well one can 
reconstruct an isolated cloud’s shape and internal structure using satellite data from 
MISR (275 m pixels), which has global coverage, rather than AirMSPI (10 to 20 m 
pixels), which is deployed only at field campaigns. We anticipate that optimization 
using this coarser data will require careful regularization. 
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Efc 8.606 10+4 km-i 

Efc 9.800 10+4 

Mean absolute relative error 4.672 10“® 

Maximum absolute relative error 8.025 10“4 


Table 3: Recovery results 



Figure 8: A simulated fly-over of the AirMSPI using SHDOM to generate the radiance 
measurements at 3 of the 9 viewing angles. The cloud extinction held is that of the 
recovered cloud. 



Figure 9: A volumetric comparison between the ground truth LES generated cloud and 
the recovered cloud. Axes are in km. 
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